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ABSTRACT 


A matrix having a high percentage of zero elements is called sparse. 
In the solution of systems of linear equations or linear least squares 
problems involving large sparse matrices significant saving of computer 
cost can be achieved by taking advantage of the sparsity. This 
Memorandum derives and describes the well known conjugate gradient 
algorithm and a set of related algorithms which are applicable to such 
problems. 

Control of accuracy is a serious problem with this class of methods. 
We plan to devote a subsequent study to methods of controlling algorithms 
of this class . 


vi 
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PART I 
Introduction 


Chapter 1 Introduction 

A matrix A is called sparse if a large proportion of its elements are zero. 
Significant savings of execution time and computer storage can be realized in 
solving systems of linear equations, Ax=b, or least squares problems, Ax — b, if 
the matrix A is sparse and if special solution methods are used which take advantage 
of the sparseness. 

Most direct solution methods such as Gaussian elimination or Householder 
orthogonal triangularization perform transformations on the given matrix A which 
significantly increase the number of nonzero elements. The two main ideas in 
developing sparse matrix methods have been 

1. Reorganize direct elimination methods to reduce the growth of the 
number of nonzero elements. 

and 

2. Use iterative methods so that the original sparse matrix A is used 
throughout the computation. 

In Reid (1971) it is pointed out that the conjugate gradient (CG) method for 
solving a system Ax=b, with A symmetric and positive definite, is well suited 
for use when A is sparse. This method shares with iterative methods the feature 
of continually using the initial matrix A but it shares with direct methods the 
property of theoretically terminating after not more than n iterations. 

More specifically the major computational cost in each iteration of the CG 
method arises from the multiplication of A times a vector. If A is an n x n 
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matrix and the proportion of nonzero elements in A is p then the multiplication of 

2 

A times an n-vector requires pn multiplications and additions if multiplication 

by zero elements of A are skipped. This algorithm theoretically reaches the 

solution vector in at most n iterations so the number of multiplications and additions 

3 

would theoretically be at most pn . 

The storage required is just that needed for the sparse matrix A plus four 

2 

n-vectors. Since A is symmetric it has only approximately pn / 2 potentially 

2 

distinct nonzero elements. Such a matrix can be stored in pn locations or even 

in less space if the nonzero elements are located in some known regular pattern. 

In comparison direct solution of this problem using an efficient stable method 

2 

such as Chclesky factorization would in general require abou* - n /2 storage 
3 . 

locations and n /6 multiplications and additions. Thus the CG method will require 
less storage than the Cholesky method if p < 1/2 and will require fewer 
multiplications and additions if p < 1/6. 

These cut-off values for p should be taken only as very rough guidelines. 

Storage management for the sparse matrix method may increase its execution 
time . 

More serious is the lack of numerical stability in the CG method. If the 
matrix A has a large condition number the intermediate vectors computed by 
the algorithm which are theoretically orthogonal or A-conjugate (defined in 

H 

Chapter 2) may not even come close to having these properties. I feel that this 
means that a reliable subroutine for the CG method must include some 
monitoring of the error generated in intermediate quantities. 

The purpose of this report is to collect in one place, and with some 
consistency of notation, the statements and theoretical justifications of the 
conjugate gradient algorithm and a number of other algorithms having very 
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similar characteristics with regard to mathematical theory, operation counts, 
and storage requirements. We subsequently plan to produce Fortran subroutines 
for some of these algorithms and study particularly the effectiveness and 
reliability of various techniques for monitoring accuracy and testing for 
termination in these subroutines. 

The CG algorithm was invented independently and simultaneously (~ 1951) by 
M, R. Hestenes and E. Stiefel [see Hestenes and Stiefel (1952)]. The paper by 
Craig (1955) which discusses methods of this type also references related work 
by Fox, Huskey, Wilkinson, Lanczos, Forsythe, and Rosser in the 1948-1952 
period. 

After providing some mathematical background in Chapter 2 the CG algorithm 

is presented in Chapter 3. In Chapters 4 and 5 algorithms are given which 

T T 

result directly from replacing the matrix A in the CG algorithm by A A or AA 

T 

respectively. The use of A A covers the case of a least squares problem 

T 

whereas the use of AA allows solution of consistent systems, Ax=b, in which A 

is not necessarily symmetric or positive definite (or even square). 

T T 

This replacement of A by A A and by AA occurs in Craig (1955) and is 
also treated in Faddeev and Faddeeva (1963). 

More recently Reid (1971) made a strong case for the usefulness of the CG 
method in large sparse problems. Interest in sparse problems also stimulated 
Paige (1972) to derive two algorithms of similar character. We will refer to 
these algorithms as PAIGE-I and PAIGE-II. Paige's derivation of these algorithms 
is based on a bidiagonalization method given by Golub and Kahan (1965) which has 
its mathematical roots in a method of Lanczos (1950) for tridiagonalizing a 
symmetric matrix. We present Paige's least squares algorithm, PAIGE-II, in 
Chapter 6, where we call it ITLS to distinguish our particular statement of the 
algorithm. In Chapter 7 we show that ITLS theoretically generates the same 
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sequence of approximate solution vectors as the algorithm of Chapter 4. The 
intermediate steps are sufficiently different however to make it of interest to 
investigate the numerical performance of each of these two algorithms. 

In Chapter 8 we present an algorithm ITC which is a subset of Paige's algorithm 
PAIGE-I. The algorithm PAIGE-I includes provision for handling a least squares 
problem. It appears to me that for least squares problems this algorithm is 
not competitive with PAIGE -II or the least squares algorithm of Chapter 4 and 
thus I have presented only the subset of PAIGE-I which handles a consistent 
system of linear equations. Paige noted that this subset of PAIGE-I is 
equivalent to the algorithms given by Faddeev and Faddeeva (1963) and by Craig 
(1955) which are described in Chapter 5 of this report. This equivalence is 
verified in Chapter 9. 

Finally in Chapter 10 we give an algorithm due to C. C. Paige and M. A. Saunders 

(personal correspondence, 1972) for the symmetric consistent problem, Ax = b 

This algorithm is called SYMMLQ by Paige and Saunders. Note that whereas 

the CG algorithm requires only one matrix-vector multiplication per iteration 

the other algorithms discussed in Chapters 4-9 each require two matrix- vector 

T 

multiplications per iteration or else the preliminary computation of A A or 
T 

AA . The algorithm, SYMMLQ, requires only one matrix-vector multiplication 
per iteration and thus is nominally the most economical method described in 
this report for the indefinite symmetric consistent problem. This algorithm is 
however notably more complicated than the other algorithms in this report. 

Saunders has also developed a modification to Paige's least squares algorithm, 
PAIGE -II, (our ' rapter 6) however we have omitted this as it is more complicated 
than Paige's algorithm and in a few test cases which we ran its sequence of 
approximate solution vectors was approximately one step behind the sequence 
generated by Paige's algorithm. 
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Chapter 2 Mathematical Background 


Two real n- vectors x and y are mutually orthogonal if their inner product , 
T T 

denoted by x y or y x, is zero. They are orthonormal if they are mutually 
orthogonal and are each of unit euclidean length, i. e. j|x|| sfx'^'x) 1 =1 and 

fly I! s (y T y ) 1/2 = l. 


A generalization of orthogonality is the notion of A-conjugacy where A is 

• T 

a symmetric matrix. Two vectors x and y are A- conjugate if x Ay (or 

T 

equivalently y Ax) is zero. This notion of A-conjugacy is most commonly 

defined only for a positive definite symmetric matrix A since then one has the 

T 

convenient property that x Ax > 0 for all x ^ 0. 

We will use the notion of A -conjugacy under the weaker assumption that A 

is nonnegative definite symmetric matrix but limit the vectors being considered 

T 

to those lying in the row space of A. For such vectors the property x Ax > 0 
for x 4 0 still holds. 

If the set of vectors Yf = {v^, . , . , v^ } are mutually orthogonal and a vector 
y( i"^ i ) ii near iy independent of the set"yC then a vector v^ + ^ orthogonal to the 
set can be defined by 


( 2 . 1 ) 


v (i+l) _ y (i+ 1 ) « v 


(l)T v (i+l) 

v^V 1 ) 



v (i)T v ( i+1) 

V^vW 



This is the formula of Gram-Schmidt orthogonal ization. 
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A similar formula exists for extending a mutually A- conjugate set of vectors. 
Thus it the set of vectors - [u^, . , . , u^} are mutually A- conjugate and a 
vector Z < iH > is linearly independent of the set then a vector u^ A-conjugate 
to is defined by 


( 2 . 2 ) „< i+1 > = 


u"> T A a <‘ +, > m 
u“> T A u<‘) 


U) 

u' , » t Au«> U 


if all of the denominators are nonzero. Nonzero denominators are assured if A 

is positive definite symmetric or if A is nonnegative definite symmetric and all 

of the vectors of the set lie in the row space of A. 

The algorithms to be described in this report have the common feature that 

in constructing sequences of mutually orthogonal (or mutually A-conjugate) vectors 

(i+1) (i+1) 

the new linearly independent vector y' ' (or z' ') will be constructed in such 
a way that it is already orthogonal (or A-conjugate) to all but one or two of the 
vectors in the set /. (or &), This permits economy in storage and in 
computation time since only the most recent one or two of the vectors in the 
set (or U^) need to be retained in storage and only the terms involving these 
one or two retained vectors need to be computed in Equations (2. 1) or (2. 2). 
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PART II 


The Conjugate Gradient Algorithm and Variations 


Chapter 3 Solving a Consistent System, Ax=b, where A is Symmetric and Nonrv'a&tivp 
Definite 

Let A be an n x n symmetric nonnegative definite matrix and let b be an n~ 
vector in the column space (range space) of A. We wish to find an n-vector x 
satisfying 

(3.1) Ax « b 

If Ranh (A) < n the solution of Problem (3. 1) is nonunique. In this case there 
is a unique solution vector, x, in the row space of A. This vector x is the 
minimal length solution vector for the problem and is the solution vector which 
the algorithm to be described constructs. 

The algorithm to be described is the conjugate gradient method due to Hestenes 
and Stiefel (1952). Presentations of this mtthod appear in Faddeev and Fadeeva 
(1963, pp. 392-405), Beckman (I960), Reid (1971), Fox (1965, pp. 208-214), 
and Householder { 1964, pp. 139-141). This method is usually described under the 
assumption that the matrix A is positive definite although, as will be seen from 
the discussion to follow, the theory of the method is also valid for a nonnegative 
definite matrix if it is assumed that b is in the column space of A, 

Assume the existence of an integer k( 1 s k s n), matrices V , and D, , 
and a k-vectcr p such that 

(3.2) V - [v (1) , . . . , v (k) ] 
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(3.3) 


D = Diag (d j , . . . , dj ] cL > 0, i=l,...,k 


(3.4) 

V T AV = D 

and 


(3.5) 

A - ■. r* 

x - Vp 


Note that Equation (3.4) implies that the vectors are mutually A-conjugate. 

If such matrices V and D are available Problem (3, 1) can be attacked as 

T 

follows: Left multiply Equation (3, 1) by V obtaining 

(3.6) V T Ax * g 
where g is defined by 

(3.7) g * V T b 
Introduce the change of variables 

(3.8) x = Vp V 

in Equation (3. 6) obtaining 

(3.9) V T AV P = g 
which due to Equation (3.4) may also be written as 
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{3. 10) 


Dp = g 


T 

Thus Problem (3. 1) could be solved by computing g - V b, solving the diagonal 
system of equations Dp = g for its solution vector p, then computing the solution 
vector x for Problem {3. 1) as x = Vp. 

The algorithm to be described constructs the A-conjugate vectors v^ one at 
a time and as each such vector is produced its contribution to g, p, x, and the 
residual vector r is determined. Thus only one of the vectors v^ needs to be 
maintained in storage at any one time. 

It will be convenient to define auxiliary vectors 

(3.11) w (i) = Av (i> i= 1 , , . . , k 


successive approximations to the solution vector 


(3. 12) 



0 


(3.13) 


,<i) = 


= S 


j=l 


(j) (i-1), (i) 

v' J 'p. = x' + v p. 


i= 1 , . . . , k 


and corresponding residual vectors 


(3.14) 



10 
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(3.15) 


r' ' = b-Ax' « b - E Av'-"p. 

J 


= b - 2 w^p. = ^ - w^'p. i=l, . . . , k 

j=l J 


It happens that the residual vectors r^, r^, . . . occurring in this algorithm 
are mutually orthogonal. The algorithm alternates between producing a vector in 
the orthogonal sequence r^, r^\ . . . and a vector in the A-conjugate sequence 
v^, v^, .... The various relations which exist between these two sequences 
permits the algorithm to be remarkably concise. 


(3.16) Algorithm CG The Conjugate Gradient Algorithm for Solving a Consistent 


System Ax=b where A is Symmetric and Nonnegative Definite 
[Due to Hestenes and Stiefel { 1952) ] 


Description 

( 0 ) _ ( 0 ) , ( 1 ) , 

x' ':- 0 , r' ':=b, v' ':=b 

If b=0 set i; = 0 and go to Step 13 


w li, :=Av (i > 

(i) (i“ 1) (i) 

r' ':=r' -w' 'p. 

Theoretical termination test: If r^=0 go to Step 13 
Practical termination test: If ||r^ij is sufficiently small 
go to Step 13 
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Step 

Description 

10 

(i+1 ) (i) , 

v' ': = r' '+v 

11 

i:=i+l 

12 

Go to Step 4 

13 

k: s i 

14 

Stop 


Figure (3. 1) is provided as an aid to understanding Algorithm CG. All 

vectors in the first column of Figure (3. 1) lie in the same one-dimensional 

} £-1 
subspace, « j . For general l > 1 if the vectors b, Ab, . . ■ , Ab are linearly 

independent let yf^ denote the X-dimensional subspace spanned by these vectors. 

Then the first k vectors in each row of Figure (3. 1) are also linearly independent 

and span the same subspace 

To verify that the algorithm CG is mathematically correct we must show 
that the denominator in Step 5 is positive for i £ k, that Step 5 defines components 
p^ satisfying Equation (3. 10), and that the vectors v^ produced at Step 10 are 
mutually A- conjugate . It will also be seen that the residual vectors are 
mutually orthogonal. 

Assume Algorithm CG has been executed for i=l, . . i-\, that the set of 

vectors {r^, ...» r^ are mutually orthogonal and the set of vectors 

tv' 1 ) v< X >] are mutually A- conjugate. With this assumption of A-conjugacy 

( i) T ( i) 

we include the assumption that v' ' Av > 0, i=l, . . . , l. 

We also assume that 
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*f. = Span{b, w^, . . . , 
s Span {r^, , . . , 

= Span {v* 1 * v^} i~l , . . . , Jt 


The reader may find it convenient to refer to Figure {3.1) and consider Z-2 

for definiteness. The quantities to the left of the circles result from earlier 

iterations and the circled quantities will be computed during the second iteration. 

Consider now the quantities computed during the iteration, The denominator 
(£)T (l) 

in Step 5 is v' ' Av' ' which by assumption is positive. 

It must be verified that p^, computed at Step 5 satisfies Equation (3. 10), i.e. 

that 

(3. 17) Pjl = g^/d^ = (v^ )T b)/(v U,T Av U) ) 

The denominator in Equation (3. 17) is clearly identical with the denominator in 
Step 5. As to the numerators we have 

j=l J 

= v U)T r U-i) 

= (r U " 1 ) +v ( ^ 1 ) ^_ 1 ) T r (^l) 

= r U-i)T r U-i) 


14 
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Steps 6 and 7 are clearly consistent with Equations (3. 13) and (3. 15). If 
= 0 the algorithm terminates, setting k=i,. Otherwise with r^ ^ 4 0 we proceed 
to verify that r^ is orthogonal to the subspace Using the basis 

{r^, . . . , r^ of it suffices to verify that r^^r^ = 0 for i=0, X,-l. 

(3.18) r U)T r (^ = r (i)T r (i-l)_ r (i)T w (^ 


Both right side terms are zero for i=0, . . . , 1-2. For i=f-l substitute the 
definition of from Step 5 into Equation (3. 18) obtaining 

U-1)T U) „ (f- 1 )T { i-l) 

= 0 


since 


U)T (j t) P (l-l) , (je-l) Q i T (&) <£-l)T 

r ' w' = [ r' + v' ® i - 1 w = r' w 


U) 


Next it must be shown that v^ +1 ^ defined at Step 10 is A-conjugate to the 
mbspace We use the basis (v^ ...» v^} for and verify that 

^> T Av^ + 1 ) = 0 for i=l, 


(3.19) 


,(i)T Av (^D =v (i)T Ar (i) +v (i)T Av (X ) 3 
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( l ) T 

For i=l, ...» £-1 both right side terms are zero; the first because v v ' A is in 

•> i ji\ [i) ? 

and thus is orthogonal to r' ' and the second because v' ' is in ®^_j and thus 

/ i\ 

is A-conjugate to v v For substitute the definition of 8^ from Step 9 into 
Equation (3.19) obtaining 


(3.20) 


^U)T Av U+1) = v U)T Ai .U) , (v U)T Av U) )(r U )T r U) 


r (je-i)T p (je-i) 


- ,„( JJ)T (&) , -1 (i) T (JL) 

- w' r + p^ r' ' r 

= Cw (W + P 2 1 r UI 1 T f a) 

= p -l r U-DT r (/) = 0 


Finally we verify that v ^ + ^^Av^^^ > 0. Let h denote the rank of A. 

Since A is nonnegative definite there exists an h x n matrix R of rank h such that 

T 

A a R i R 


and clearly the row space of R is identical with the row (and column) space of A, 
Let Ct denote the row (and column) space of A. From Steps 1, 7, and 10 all of 
the vectors r^ and v^ produced by the algorithm lie in the subspace & . Since 
v U+1) is constructed at SteplOas the sum of the nonzero vector r^ and a vector 
v^8. orthogonal to r^ it follows that v^ + ^Vo. Since v^ + ^ is a nonzero vector 
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in the subspace d the vector Rv^ ' 1 ^ ^ must also be nonzero. It follows that 
[Rv< 1+1 >] T Rvl 4+1 '>0o r equivalently > 0, as was to be jhown. 

As is apparent from the different expressions derived for p^ and for B. 
there are a number of different ways in which the CG algorithm could be 
implemented. There are also different trade-offs possible between storage used 
and counts of arithmetic operations. Reid (1971) discusses over a dozen such 
variations. The form in which we have stated Algorithm CG is the one 
preferred by Reid. 

Theoretical termination of the CG algorithm occurs when r^=0 with r^^O 
for i=0, . . . ,k-l. Referring to Figure (3. 1) we see that this means that the 
subspaces . . . , 9 /^, are all different but that = 9^. Equivalently this 

means that the vector A^b lies in the subspace ^ spanned by {b, Ab, . . . , A^ *b}. 

Tiiis implies that b is representable as a linear combination of some set 

of k eigenvectors, say (£<*> £< k >], of the nonnegative definite symmetric 

matrix A and is not representable as a linear combination of any smaller set 
of eigenvectors of A. Furthermore the eigenvalues . . . , associated with 

these eigenvectors are all positive. Thus there exist nonzero numbers c. such 


k , .v 

b = E f ^ c . 
i=l 1 


and the minimal length solution vector x is representable as 


k ... 

x = E f (1, (c./\.) 
i=l 1 1 
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Note also that Span{f 


, f (k *3 - v/^ and that x lies in vJ 1 ^ but not in xf^_ j . 

Most commonly the value of k will be n. However k will be less than n if 
and only if b is orthogonal to some eigenvectors of A. In particular b will 
necessarily be orthogonal to some eigenvectors of A if A has any multiple 
eigenvalues. 


18 
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Chapter 4 Solving a Least Squares Problem Ax — b 


Let A be an m x n matrix and let b be an m-vector, We wish to find an 
n-vector x which minimizes llb-Ax^. We denote this least squares problem by 
the notation 

(4.1) Ax s: b 

We permit either m ^ n or m < n. If Rank (A) < n the algorithm to be described 
constructs the unique minimal length solution vector, £ . This solution vector is 
characterized ss being the only solution vector lying in the row space of A. 

A vector x minimizes |tb-Axjl if and only if it satisfies the "normal equations" 

(4.2) A T Ax - A T b 

T 

Normal equations are always consistent since the right side vector, A b, lies 

T 

in the row space of A which is also the column space of the matrix A A. The 

T 

ranks of A A and A are equal and their row spaces are the same. Thus if 
T 

A A is singular the unique solution of Problem (4.2) lying in the row space of 
T 

A A is also the unique minimal length solution of Problem (4. 1). 

T 

Since A A is nonnegative definite and Problem (4. 2) is consistent the 
conjugate gradient algorithm CG (3. 16) is directly applicable to Problem (4.2). 
Denote the residual vector for Problem (4.2) by 

h = A T b-A T Ax = A T (b-Ax) 
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and introduce bars on various other symbols in Algorithm CG to distinguish 
the application of the algorithm to Problem (4.2). Then Algorithm CG adapted 
to Problem (4. 1) can be written briefly as 

(4.3) Algorithm CGLS Conjugate Gradient Algorithm for the Least Squares 

(4.4) 

Problem Ax 2 * b 

*"»•. = 0, h* 0 '-. =A T b, v<‘> ; =A T b 

(4.5) 

If h^=0 terminate. 

Do for i: = l,2, .,,, until h^=0 

- A T Av^ 

(4.6) 

p.: - |lh (i ‘ 1 >|| 2 /(v {i)T w (i )) 

(4.7) 

r i 

(4.8) 

r i 

(4.9) 

0.: = |ltf l) || 2 /ltf i " 1) l |2 

(4.10) 

v* 1+1 *: = h*^ + v^eT. 

l 


One may wish to have the residual vector of the least squares problem (4. 1) 


(4.11) rW = b-A5EW 

available at each iteration, possibly for use in a supplementary termination test. 
This may be accomplished by the following revision of the algorithm fFaddeev and 
Faddeeva ( 1963, p, 403]. 
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(4. IE) Algorithm [Alternate Form of CGLSl 


2<°W. -?<°W, h‘°>:=A T b, Vl‘»: 
If h^=0 terminate. 

Do for i: = l , 2, . . . , until 17^=0 

u^: = Av (i) 

P i : = i|H (i - 1 ) H 2 /|u < i )|| 2 


= A T b 


x (i) : =x (i_1 > + v (i) p. 


— ( i) . _ _{i- 1 ) — (i)~ 

r : - r - u o. 

*i 

K< S >: = A T r<‘> 


r,i - iih (i, ii 2 /iih (i - 1 ) p 2 


This latter form of the algorithm requires more storage since one must 
maintain the current values of the two m"vectors u^ andlr^. 
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Chapter 5 Solving a Consistent System Ax-b 


Let A be an m x n matrix and let b be an m-vector contained in the column 
space (range space) of A. Consider the problem of finding an n-vector x satisfying 

(5.1) Ax = b 

We will permit either m £ n or m > n. If Rank(A) < n the solution to this 
problem is nonunique. In this case there is a unique solution vector x lying in the 
row space of A. This vector x is the unique minimal length solution vector for 
Problem (5. 1). The algorithm to be described constructs this solution vector x. 

Since we seek a solution vector x in the row space of A the solution vector x 
will be representable in the form 

(5.2) ' x - A T y 


for some (not necessarily unique) m-vector y. Making the change of variables 
T 

x = A y in Problem (5. 1) we obtain the problem 

(5.3) A^> T y = b 

This is a consistent problem with a symmetric nonnegative definite matrix 
T 

AA . Thus the conjugate gradient algorithm CG (3. 16) can be applied to solve 
Problem (5. 3). 

The resulting algorithm for solving a consistent system Ax=b may be written 

as follows where tie notation of Algorithm CG (3, 16) is changed by writing p, 8, 

_ _ _ 1 * 

w, r, q f y, and AA in place of p, 8, w, r, v, x, and A respectively, Furthermore, 


22 


JPL Technical Memorandum 33-627 


motivated by Equation (5.2) we introduce the sequence of approximate solution 
vectors 


(5.4) 


x^ = A T y^ 


(5.5) Algorithm 

y(°> ; = 0, *<°>: = 0, T«» : - b, ,<»•. 
If - 0 terminate. 

Do for i:=l, 2, . . . , until r^ = 0 

w^: = 

p.: = llr* 1 " 1 ^ 2 /^^ T w^) 

yW =? (i_I) + qWp. 


S b 


r ( 1 ), _ r^^-vS^p. 


F.: « l!r (i) || 2 /||r (i " 1 ) |l 2 

^(i+i). _ ^(1) + qU)g- 


Eliminating the intermediate vectors and y ^ we obtain the algorithm in 
the form given by Craig (1955), p. 72, except that Craig used the opposite choice 
of the sign of the residual vector. 
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(5,6) Algorithm [ Craig { 1955) 3 


3c<°>, = 0, - h. q(') : =b 

I£ -J-(O) s q t erm i na £ e> 


Do for i:=l, 2 until = 0 


P.: = ||r< i “ 1 >|! 2 /!lA T T T < i) l | 2 
x^: = *( 1_1 ) + A T q^p. 

7^: = b-Ax^^s 7^ 1 ^ - AA T q^^pjj 
F.: = ||r (i) ll 2 /|M i ' 1) ll 2 

•q;( 1+1 >: = + q^j-T. 


As is noted in Faddeev and Faddeeva (1963, pp, 403-405) this algorithm can 
be further simplified by introducing the substitution 

-(i) A T-(i) 
v' = A q' 

Note that the vectors {v^, . . . , v^J are mutually orthogonal since the vectors 
{q^i . . . ,q^)] are mutually (AA T )-conjugate, 

We call the resulting algorithm CGC. 


(5. 7) Algorithm CGC Conjugate Gradient Algorithm for the Consistent System Ax=b 
3E< 0 >:=0, r<0) :=b> *D. = A T b 


If r^ = 0 terminate 

Do for i:=l , 2, . . . , until r^ = 0 
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P t ! » l|7 (i ' 1| l| 2 /l!v (i >l| 2 
X<‘> : = 

r (i >: = r< i -‘>-Av< i >p. 

1 1 

t>V = Il7 (i) |l 2 /ll7 (i ‘ 1) ll 2 

-( 1 - 1 ) ; „- A T F (i) + -U) B 
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PART III 


OTHER ALGORITHMS 

Chapter 6 Solving a Least Squares Problem, Ax a; b 

Let A be an m x n matrix and let b be an m- vector, We wish to find an n- vector 
x which minimizes {|b-Axjl. We denote this least squares problem by the notation 

(6.1) Ax — b 

Generally one would have m £ n and Rank(A) = n in such a problem. We will 
not make these assumptions however as they are not necessary to the mathematical 
development of the algorithm to be described. 

If Rank (A) < n the solution to Problem ( 6. 1 ) is not unique. In this case 
however the problem has a unique solution vector of least euclidean length. It 
can easily be verified that this unique minimal length solution vector lies in 
the row space of A and in fact is the only solution vector for Problem (6.1) lying 
in the row space of A. The algorithm to be described constructs the solution 
vector in the row space of A and thus finds the minimal length solution vector if 
Rank (A) < n. 

Let x be the unique minimal length solution vector for Problem ( 6. 1 ). 

Define the residual vector 


(6.2) 

A v . A 

r = b-Ax 

and 

(6.3) 

b = Ax = 


Thus b can be written as 

(6.4) b = b + r 
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where b lies in the column space of A and r is orthogonal to the column space of A. 
The vector b is the orthogonal projection of b into the column space of A and will be 
referred to as the projected right-side vector . 

Suppose there exists an integer k (1 i k i min Cm, nl) and matrices 


(6.5) 


U . = [a* 1 *. . . .,u {k) ] 
m x k 


( 6 . 6 ) 


W=tv<‘> <“>] 


(6.7) 


l kxk 


"1 B 2 


«2 ' B k 

* . 


o', 


(all a. > 0 and EL > 0) 


and a k-vector, p, such that 


( 6 . 8 ) 

(6.9) 

( 6 . 10 ) 

( 6 . 11 ) 

and 

( 6 . 12 ) 


U T U = I k 

v T v - I k 

AV = UR 
A T U = VR T 


A 

X 


Vp 
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Since R is nonsingular Equations (6.9) and (6.11) imply that the column vectors 
of V form an orthonormal basis for a subspace, [/, of the row space of A. Similarly 
Equations (6.8 ) and (6. 10) imply that the column vectors of U form an orthonormal 
basis for a subspace of the column space (range space) of A. Equation ( 6. 12) 
shows that the subspace V* contains the solution vector x. From Equations (6.3 ), 

(6. 10), and (6. 12) it follows that the subspace ll contains the projected right-side 

A 

vector b. 

Assuming the availability of the matrices U, V, and R, Problem ( 6. 1 ) can be 
approached as follows: Introduce the change of variables 

(6.13) x — Vp 

in Problem ( 6. 1 ) and use Equation (6. 10) obtaining the equivalent least squares 
problem 


(6.14) 


URp ^ b 


Let U m x (m-k) a ma * :r * x which when adjoined to U forms an m x m orthogonal 
matrix: 


(6.15) 


[U:U] T [U:U] = I m 


*T 

Left multiply Equation (6. 14) by [U:U] obtaining the equivalent least squares 
problem. 


(6.16) 



P - g s 




k 

m-k 
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where 


(6.17) 


6 nn 

= [U:U] i b = 


U T b 


The least squares solution vector p for Problem (6. 16) may be obtained as the 
solution of the upper bidiagonal nonsingular system of equations 


(6. 18) 


Rp = g 


Obtaining p from Equation (6. 18) the solution vector x can then be computed from 
Equation (6. 13). 

In analogy with the conjugate gradient algorithm we wish to cast this procedure 
into a sequential form so that a succession of approximate solution vectors x^ 
and associated residual vectors 


(6.19) 


r< x) = b-Ax* 1 ' 


/ * \ i * \ 

can be computed as the successive vectors u* ' and v' ' are computed. This 
obviates the need to store old u^ and v * vectors. 

Equation (6. 18) is not directly suitable for such a (forward) sequential 
procedure since the last (lower right) element of R must be determined before 
any components of p can be computed. However using Equations (6. 13) and 
(6. 18) we can write 


( 6 . 20 ) 


x = Vp = VR^Rp = (VR'^g = Wg 


JPL Technical Memorandum 33-627 


where the n x k matrix 


( 6 . 21 ) 


W 



w<“>] 


satisfies the linear matrix equation 
(6.22) R T W T = V T 


or 


(6.23) 


r 

*1 

0 " 


r (ijTi 

w 


V (')T 

8 2 

» 

«2 

# 

4 

• 

• 

• 

* 

• 

« 

w WT 

= 

» 

t 

v (k)T 

0 

* 9 k “k 






From this matrix equation one obtains the following expressions for sequential 
computatiali of the vectors w^. 

(6. 24) w^ = 

(6.25) w^ = (v^-g.w^ ^)/a. i=2, , . . , k 


L.et denote the i th component of the m-vector g. 
the i component of the k-vector g of Equation (6. 18). 


Then for i £ k, is also 
Define the m-vectors 
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C g 2 » • * • » g^ 


o,,..,o i T 

m-i 


i-0, 


m 


and the k-vectors 


Mi) 

g 


= Cgj. 


®i’ 


0, 


0]' 


i=0, 


k-i 


k 


Motivated by Equation (6.20) define a sequence of approximate solution 

* (i) u 

vectors x by 


( 6 . 26 ) 


( 0 ) 


0 


(6.27) 


= Wg (l) = E 


j = l 


(i- 1 ) (i) 

= x' ' + w' 'g^ 


i= 1 k 


( i } 

The associated residual vectors, r' \ are defined as 

(6.28a) r^°) = b 

and 

(6.28b) 


,<« = b-Ax^ 


= b-AVR g l ' 


b-URR -1 g (l) 


[Using Equation (6. 10)] 


b-uS<‘> 
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b- i u^g. 

j=i 1 




i— 1 lc 


We next consider formulas for computing 'g and HA^r^jj, i-1, . . , , k, which 
depend upon a particular choice of v^ . The choice of is somewhat arbitrary 
as long as it is chosen to lie in the row space of A. We follow Paige (1972) in 
defining 3j and by 


(6.29) 


v ( 1 *9 2 = A T b 


where 8^ = ||A *^b 11 so that jjv^jf = 1. 
From Equation (6. 17) we have 


(6. 30) 


r r 

g = u A b 


Eeft multiplying this equation by R A and using Equation (6. 10) and (6.29) gives 


(6. 31) 


T T T T 

Rg = R A U b = V A A A b 




where e^ denotes the first column of the kxk identity matrix. 

T~ 

Writing the equation R g = 9jej in terms of its components gives 
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(6. 32) 


^2 

0 


on 




from which the components can be computed as 

(6. 33) g 1 = 6 1 /a' 1 

(6.34) g i = -(0 i /c r i )gi_ A i=2, . , , ,k 


Note that it would also be possible to compute the quantities g^. sequentially as 

g^ = u^^b (see Equation (6. 30)) but Paige (1972) reports that Equations (6. 33)-{6. 34) 

were found to give better numerical accuracy in test cases. 

In a least squares problem one does not generally expect the final residual 

vector, r = b-Ax, to be zero. The residual vector at the solution is characterized 

however by the property that it is orthogonal to all of the column vectors of A. 

T 

Thus the vector h - A (b-Ax) is zero if and only if x is a solution vector for the least 

squares problem Ax — b. It is also true that h is the negative gradient vector with 

1 . 2 

respect to x of the function ^Hb-Axl 1 . Define 


(6.35) h( i )=A T r< i )=A T (b-Ax< i )) 

= A T (b-ui (i) ) 

= v (i) 0 1 -VR T ^ i) 
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= V(e (1) 8 1 -R T g (l) ) 
’ • VeU+1)s i + l‘!i 


i=l , . . . , k-1 


while for i=0 and i=k we have h^=A*^b-v^ ^0 ^ and h^-0. It is of interest to note 
that the vectors h^, . . , , h^ k ^ are mutually orthogonal. 

The quantities which may be useful in monitoring the progress of the 

algorithm are thus expressible as 


(6. 36) 


Y 0 “ 8 1 


(6.37) 

(6.38) 


Y i - i p i+iBj I i=1 » .... k-i 

Y k = 0 


In practice k is generally not known in advance and will in fact be defined as the 


first value of i. for which 8^ " 0. 


(i) 


We turn now to the determination of the quantities u' ' and for i=l, . . . , k and 


,<‘1 

(6.39) 


v' 1 * and B^ for i=2, . . , , k. From Equation (6. 10) we obtain the equations 

u^-Av* 1 * 


(6.40) 


u (l) a. = Av (l) -u (l_1) B. i=2, . . . , k 


and from Equation (6. 11) the equations 


(6.41) 

(6.42) 


v^p. = A T u (i " 1) -v (i “ 1) a._ 1 i=2, . . .,k 

0 * A T u< k > - v< k >a, 


If B j and v^ are defined by Equation (6. 29) and the scalar quantities B. and a?, 
are determined so that the vectors v^ and u'^ respectively have unit euclidean 


length (as is required by Equations ( 6. 8 ) and ( 6. 9 )) then Equations (6. 39)-{6. 41) 
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determine all of the remaining vectors , i-2, . . , , k, and i«l, . , . , k. 


Collecting these various formulas together leads to the following algorithm. 


(6.43 ) Algorithm ITL-S for iterative solution of the least squares problem 
Ax £: b. [Due to C. C. Paige(l972) pp. 21-22.3 


Step No. 

1 

2 

3 

4 

5 

6 
7 


8 

9 

10 

11 


12 


13 

14 

15 

16 

17 

18 


Description 


x<°U 

go -" 1 

i; -l 


A T b if i=l 

A T u (i-l)„ v (i-l) a ^ 

l- 1 


B i :=[|v <i> ll 

v i- l t= I ®i ®i- 1 1 

Theoretical termination test: If 0. = O go to Step 17. 
Practical termination test: If either 8^ or j is 
sufficiently small go to Step 17. 

v< 1 >:=r< i >/B. 

1 
r 

Av^ 1 ) if i=l 
[av^^-u^ ^0. if i > 1 


Mi) 

u w< = 




w 


(i). _ 


v^Vnfj if i=l 

|( v ^ 1 )- w ( 1-1 )j5 i )/ ai if i > i 


(i) (i-l) x (i) 
x' ':=x' '+w' 'g^ 

i: ~i+ 1 

Go to Step 4 

k:=i-l 

Stop 
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It must be verified that the vectors v^ and produced by this algorithm 
have the orthogonality properties specified by Equations { 6. 8 ) and (6.9 ) and 
that all of the numbers o\ defined by this algorithm are positive. 

Assume that i~\ iterations of Algorithm (6.43) have been executed producing 
positive numbers 8^ and a^\ i— 1 , . . . , i-l, aset of mutually orthogonal unit 
n-vectors {v^,.,,,v^ ^}, and a set of mutually orthogonal unit m-vectors 
{u^, . . . , u^ ^ } . It is obvious from Steps 4 and 8 that all of the vectors v^ 
lie in the row space of A and from Steps 9 and 11 that all of the vectors u^ lie in 
the column space (range space) of A. 

Consider the quantities computed during the X*"* 1 iteration. 

If 8 ,, computed at Step 5 is zero then tne (theoretical) algorithm terminates, 

At 

setting k=j£-l. In this case we have and which means (see Equation 

(k) 

(6. 35)) that the most recently computed approximate solution x' ’ was in fact the 

unique minimum length solution £ for the least squares problem, Ax — b. 

If B „ computed at Step 5 is not zero, i. e. , B . > °’ then we raust verif y that 

^ /v/ o\ m 

the vector v' ' previously computed at Step 4 is orthogonal to v' , i=l, . . . , X-l. 

Using the formula of Step 4 with i > 1 the inner products to be investigated are 

(6.44) v (i)T 7<*> „ t (i|T A I u (H). v (i)T l ((-l)^ | 

r (0 , n T (i-1) (i)T (i- 1) 

= [u' c*j+u' 'B^J u' -v' v' a i - 1 

(i)T (£-1) , (i-l)T (X-l) a (i)T (f-1) 

= u u Of.+u u B.-v v ot . , 
l i l-l 

For i < l- 1 each of the three right side terms in Equation (6.44) is zero 
because of the assumed mutual orthogonality of {u^\ . . . ,u^^3 and the assumed 
mutual orthogonality of [v^ v^ ^}, 


36 


JPL. Technical Memorandum 33-627 



For i=£-l Equation. (6.44) becomes 


v {je-i)T~U) 


v i + °- 


a 


JL~ i 


= 0 


'*'(£) (i) 

which completes the verification, that v' is orthogonal to v' , i~2, 1. The 

( 1 )T ( £,- 1 ) 

verification that v' ‘ v' -0 is equally straightforward. 

Similarly it can be shown, that computed at Step 9 satisfies u^ T u^ = 0 
for i=l, . . . , 1 . We further assert that u' ’ is not zero and thus that a ^ > 0. 

Assume the contrary. Then 


= Av {£,) -lAv {i ~ l) -u {S, ~ 2) & ^^6 

4 ... 

= ■ • • = E c.Av' 1 ' 
i=l 1 

4 ... 

= A E c . v 
i=i 1 


where the coefficients c^, i=l, . ..,4 are nonzero. Since the vectors v^, 

i=l, . . . , 4 constitute an orthonormal basis for a subspace of the row space of 
4 ,.v 

A the vector z = E c.v' . must be a nonzero vector lying in the row space of A. 
i=l 1 

Such a vector must satisfy Az^O contradicting Equation (6.45). We conclude 
that u^^Q and a ^ > 0. 

This completes the verification of the theoretical algorithm (6.43). in practice 
there remains the problem of fully specifying a satisfactory termination test 
at Step 7. 
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Some estimate, say e., of the norm of the round-off error vector associated 

i) 

with the computed vector v' could be computed. Then the number 3. could be 
regarded as being sufficiently small for termination iffl. S e,, 

Since represents an evaluation of A^"(b~Ax^ one might define 


(6.46) 


«,= llAlKUbll + llAlMIx 1 


(ill 


and terminate the algorithm when j ^ T|uh ^ where 1] denotes the relative 
machine precision, (Define T) to be the smallest number such that the computed 
value of 1 + T] is distinguishable from 1.) 

More complex algorithmic logic might be 'needed. Thus if Fb £ but 
Y^_ j» T|cih _ j then rather than accepting x^ ^ as the solution it might be useful 
to restart the algorithm, attacking the modified problem 


(6.47) 


Adx h= b-Ax 


(i"l) 


f « | i 

to obtain a correction vector dx to be added to x' . 

We intend to study the problem of termination tests for this algorithm and 
treat the subject in a subsequent report. 
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Chapter 7 The Theoretical Equivalence of Algorithms CGLS and ITLS 


We will show that the sequence of approximate solution vectors x^ generated 
by Algorithm CGLS (4.3) is identical (theoretically) with the sequence of approximate 
solution vectors x^ generated by Algorithm ITLS (6. 43), 

It will be convenient to state the relationships represented by Algorithm 
CGLS in matrix form. In terms of the quantities defined in Algorithm CGLS define 
the matrices 


(7. 1) 

T> 

7> 

i i 

II 

(7.2) 

H =[h (0) , ...,h (k ‘‘h 

(7.3) 

C = Diagt||Av (1) ||, ...,|]Av (k) ||J 

and 

(7.4) 

= Diag [dj , . . . ,d k ) 

F= Diag{||h< 0, || ilh< k " l) !l3 

— Diag [fp* ...» fj^_j } 

-*•/ O 

Then from the orthogonality of the vectors h' we have 

(7.5) 

— T — 2 

H i H = F 

T 

and from the (A A)- 

— (i) 

•conjugacy of the vectors v' ' we have 

(7.6) 

V T A T AV = D 2 
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Following Householder ({1964) pp. 2 and 139-141) define 

i 


{ (7.7) 

J 


Then Equations (4.8) and (4. 10) can be written respectively as 

I 

-i 

] 

\ (7.8) H(I-J) =A T AVF 2 D" 2 

j 

| and 

’if 

f 

(7.9) VF' 2 (I-J T )F 2 = H 


In terms of these quantities we now define quantities, distinguished by a 
caret, which will be shown to satisfy the relationships of Algorithm ITLS. 

A _ - l 

V = HF 

A. — - 1 

U = AVD 

R = DF (I-J 1 )F 

A - 1 

W = YD 


Note that R is upper bidiagonal. 
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Using the definitions (7. 10) - (7. 12) and the Equations (7, 5) - (7. 9) it can be 
directly verified that the equations 


(7. 14) 


{7. 15) 


(7. 16) 
and 

(7. 17) 


9 T $ = I 


A 'f A 
U A U = I 


A A A 

AV = UR 


Trt A AT 

A a U = VR 1 


are satisfied. 

A 

From Equation (7. 10) the first column vector of the matrix V is equal to 

h^/jjh^P or equivalently A ^ b/ IfA^b l|. This condition along with Equations 

A A A 

(7. 16) - (7. 17) assure that the matrices U, V, and R are identical with the 
matrices U, V, and R which would be generated by Algorithm ITLS in solving the 
least squares problem Ax — b. 

Furthermore by use of Equation (7. 9) it can be directly verified that the 
matrix $ of Equation (7. 13) satisfies 


(7.18) 


wfe = 0 


In accordance with Equation (6. 30) define 


(7.19) 


A *T, 
g = U b 
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Let denote the i th column vector of W and let denote the i component 
of g. We wish to show that the correction which is added to the approximate 

solution vector at Equation (4. 7) of Algorithm CGLS is identical with the correction 

^ j » \ 

vector w' which is added to the approximate solution vector at Step 14 of 
Algorithm ITLS. Thus we wish to show that 


(7. 20) 


-U)- A(i)A 
v' 'p^ = w' 'g. 


i= 1 , . . . , k 


From Equation (7. 13) we have w^ = v^/d. and from Equations (4. 6). (7. 3), and 
2 2 

(7 . 4) we have p^ = f . _ j /d^ . Thus Equation (7. 20) may be established by proving 
that 


(7. 30) 


f i = d.g. 

1-1 i b i 


i“ 1 1 • . • j k 


Introduce the k~ vector e = [ 1, . . , , 1 


] T 


so that Equation (7. 30) can be written as 


{7.31) 

F e = Dg 



= DU T b 

[ Using Equation (7. 19)1 


= V T A T b 

[Using Equation {7. 11)] 


= V T h<°> 

[Using Equation (4.4)] 


T 

Left multiply Equation {7.8; by V and use Equation {7. 6) obtaining 
(7.32) V T H(I-J) = F 2 
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2 

Substitute this expression for F into Equation (7,31) obtaining 


(7.33) V T H(I-J)e = V T h* 0) 

as the equation to be verified. This equation is clearly true since 
(I-J)e = C 1 . 0, . . .,0] T . 

This completes the verification that the algorithms CGLS and ITLS produce 
the same sequence of approximate solutions. It follows that the vector h^ 

= A^[b-Ax^"| of Equation (4.8) is identical with the vector h^ = A^[b-Ax^] 
of Equation (6. 35). 


JPJL Technical Memorandum 33-627 


43 



Chapter 8 Solving a Consistent System Ax=b 


Let A be an m x n matrix and let b be an m-vector contained in the column 
space (range space) of A, Consider the problem of finding an n-vector x 
satisfying 

(8. 1) Ax - b 


The most common case of a consistent system of linear equations would be 
the case in which the matrix A is square and nonsingular. Also of practical 
interest is the case of a full rank underdetermined problem, i. e. , m < n and 
Rank (A) = m. The algorithm to be described permits either m ^ n or m > n and 
does not require any restriction on the rank of A. 

If Rank(A) < n the solution to Problem (8.1 ) is nonunique. In this case there 
is a unique solution vector x of minimum euclidean length. This minimum length 
solution vector is characterized by being the only solution vector for Problem (8.1 ) 
lying in the row space of A. The algorithm to be described constructs this 
solution vector, x. 

Assume there exists an integer k (1 s: k £ min{m, n}) and matrices 

< 8 - 2 > = u<k, 3 

< 8 - 3 > V n X k=Cv (1 > v (k) ] 

(all Oj > 0 and 8^ > 0) 


(8.4) 


■'k x. k 


a. 


8o of' 


8 k °k 
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and a k-vector, p, such that 


(8.5) U T U = I k 

( 8 . 6 ) V T V = l k 

(8. 7) AV = UL» 

(8.8) A T U = VL T 

and 

(8.9) x = Vp* 

Since L. is nonsingular Equations (8. 6 ) and (8.8) imply that the column 
vectors of V form an orthonormal basis for a subspace Y' of the row-space of A. 
Similarly Equations (8. 5 ) and ( 8. 7 ) imply that the column vectors of U form an 
orthonormal basis for a subspace "U . of the column space (range space) of A. 
Equation (8.9) shows that the subspace Y" contains the solution vector x. From 
Equations (8.1 ), (8.7 ), and (8.9 ) it follows that the subspace ~~U contains the 
right-side vector b. 

Assuming the availability of the matrices U, V, ond JL., Problem ( 8. 1 ) can be 
approached as follows: 

Introduce the change of variables 

(8.10) x = Vp 

in Problem ( 8 . 1 ) and use Equation ( 8 . 7 ) obtaining the equivalent problem 


JPL Technical Memorandum 33-627 


45 



( 8 . 11 ) 


ULp = b 


T 

Left multiplying this equation by U and using Equation (8.5) gives the k x k 
nonsingular lower bidiagonal system 

( 8. 12) Lp = g 

where g is the k-vector defined by 

(8. 13) g = U T b 

A computational algorithm can thus be based on the computation of g using 
Equation (8. 13), solving for p in Equation (8. 12), h.n.1 finally computing x using 
Equation (8.9 )■ These steps are all directly amendable to being organized in 
a sequential form which uses the vectors u^ and v^ as they are produced. 

Assuming the nontrivial case of b ^ 0 we follow Paige ( 1972) in defining 
u< J > by the equations 

(8. 14) 

(8. 15) 

Thus u* 1 * is in the column space (range space) of A since it was assumed that 
this was true of b. 

The other vectors and v^ and the elements o\ and EL of the matrix L are 
sequentially determined by the unit length requirements of Equations (8.5) 
and ( 8.6 ) and the following equations which follow directly from Equations 
( 8.7 ) and (8.8 ) 


= IN! 

u (1) = b/Sj 
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(8. 16) 

u^0. - Av^ 1 * 1 ^u^’ i ^ i _j 

i"2f • « i f k 

(8.17) 

0 = Av^ k U* k ^ k 


(8. 18) 

(1) _ a t (1) 

v = A u 


(8.19) 

(i) .T (i) (i“ 1) 0 

v' = A u' '-v 'Bj 

i-2, . . . , k 


Comparing Equations (8.16) and (8. 17) we note that the integer k is generally 
not known a priori but is determined as the first value of i for which Av^^-u^^a, = 0. 

With u<*> chosen in the column space of A it follows from Equations (8. 1 6) > 

(8. 18), and (8- 19) that all u^ will be in the column space of A and all will be 
in the row space of A. It will subsequently be verified that vectors u^ and v^ 
produced in this way necessarily satisfy the orthogonality conditions of Equations 
(8.5 ) and ( 8.6 ). 

With u<*> defined by Equation (8.15) the vector g defined by Equation (8. 13) 
is representable as 

(8.20) g=8 1 e* 1) 

where e^ denotes the first column vector of the k x k identity matrix. Using 
this expression for g Equation (8. 12) becomes 

(8.21) Lp = 8je^ 

which permits the components, p., of the solution vector $ to be expressed as 
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( 8 . 21 ) 


pj = 9 j /o^j 


(8.22) Pi s “(0 i /a' i )p i _ j i"2, . . . , k 

Define the sequence of k-vectors 

(8.23) p <l} = C Pl> ... tP ., 0, ... ,0 ] T i=0, , . . , k 

k-i 

Define the sequence of approximate solution vectors 


(8.24) 

*<°> = 

0 


(8.25) 


Vp' 1 ) 

= 2 v^p. = x* 1_1 * + v^p. i=l, . , . , k 

j = l 1 

Associated 

residual vectors 

are defined by 

(8.26) 



r^ = b ~ Ax^ i=0, . . . , k 

and may be 

expressed 

as 


(8.27) 


r<°> 

= b = u(‘> Bl 

(8.28) 


r«> 

= b>AVp (l) = b-UL P (l) 
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~ b-U 


0 

R 

0 

♦ 

■ 

0 


i+l p i 


♦- row i+1 


= b 


{1) .(i+l) n 

u B i u B i+i p i 


(i+i) 

-u 


8 i+l p i 


i=l, . . . » k-1 


(8.29) 



If we introduce the additional definitions, pg - -1, and the norm 

of the residual vector can be expressed as 


(8. 30) 




i=0, 


k 


These considerations may be organized into a computational algorithm as 
follows: 
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(8.31) Algorithm ITC For Iterative Solution of Consistent Linear Equations 

[Given by Paige (1972) pp. 21-22] 


Step Number 
1 
2 

3 

4 

5 

6 

7 

8 
9 

t 

10 

11 

12 

13 

14 

15 

16 
17 


Description 

* (0 >:=0 

p 0 : =' 1 
i: — 1 


b if i=l 

. (i-1) (i-1) ...... 

Av '-u “'i-i if 


~(i) 

u' = 

B i : = |IS (i, ll 

p i-l : = l B i p i-l I 

Theoretical termination test: If 8^=0 go to Step 16. 

Practical termination test: If either 8. or p. , is 

i 1-1 

sufficiently small go to Step 16. 

A T u^ if i=l 
A T u^-v* 1-1 *fl. if i > 1 


Mi) - 

v w : = 


v<‘>: . 

*1 

i: = i+l 

Go to Step 4 
k:=i- 1 
Stop 


It must be verified that all of the vectors and produced by this 
algorithm have the orthogonality properties specified by Equations (8.5 ) and 
(8. 6 ) and that all of the numbers defined by this algorithm are positive. 
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Assume these conditions are satisfied for i=l, . . . , £-1. Consider the 


th. 

quantities computed during the £ l iteration. 

If 0^ computed at Step 5 is zero the algorithm terminates setting k = l- 1. 

Thus and P k =0 which (see Equation (8. 30)) implies that the current 

(k) 

approximate solution vector x' ' is actually the unique minimal length solution 

A 

vector x. 

If 8^0 then > 0 and the orthogonality of relative to u ^ \ . . . , u^~ ^ 
must be verified 


(8. 32) 


(i)T~U) - (i)T . (£-1) (i)T (£-1) 
u' u' - u w Av' '-u w u' a s,-l 


r (i) , (i-1). -,T (£-1) (i)T (£- 1) 

= [v a.+v' '9.J v' '-u v u a, , 

1 1 


i=2 /-I 


Fori < £-1 each term of this final expression vanishes while for i = £-1 the final 

/ 1 \ * / A \ 

expression reduces to - 1 = The verification that u' ' u' -0 is 

similarly straightforward. 

n 

Similarly it can be verified that the vector v' ' defined at Step 9 satisfies 
v (i)T^;(.f)_Q £ or It remains to be shown that ot ^ defined at Step 10 

is positive. 

Suppose a.=0. Then 
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(8.33) 


n „ '-'U) - A T ( j£) (j J-1V. 

0 = v' ' = A u v' *15 




« £ c.A T u< 1 > = A T S c.u'‘) 


i=l 1 


i=l 1 


Here the coefficients c. are all nonzero and the vectors u^ f . . . , u^ constitute 

an orthonormal basis for a subspace of the column space (range space) of A. 

£ (i) . 

Thus the vector z = E c.u' ' is a nonzero vector in the column space of A. It 
T i=1 1 

follows that A ziO which contradicts Equation (8. 33} We conclude that a ^ > 0. 

The practical termination test at Step 7 might be implemented as a comparison 
of (3. with some computed estimation of the norm of the round-off error vector 
associated with the computed vector u^ or a comparison of ^ with T|( |jb || + |!,A|j * |ix^ 
where T) denotes the relative machine precision. 


* 

U 
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Chapter 9 The Theoretical Equivalence of Algorithms CGC and ITC 
Using the notation of Algorithm CGC (5.7) define 


(9. 1) 

v (i) = (-i) i " 1 v ,i) /||V (i) || 

i-1 » • 

. • » k 

(9.2) 

u (i > . (-l) i ' 1 r (i ‘ 1 >/llr (i ' 1) !l 

i=l, . 

. . » k 

(9.3) 


i=l, . 

. . , k 

(9.4) 

9, = llbll = l|7 (0) |l 



(9.5) 

IF w - 1, l!-ltr ,i - 1) l!/lt? (i ‘ 2) ll 2 




r< 

H 

if 

•H 
• H 

II 

i=2, . 

. . , k 

(9.6) 

P 0 = “I 



(9.7) 

P i =(“i) i ’ 1 |l7 (i '- 1) |l 2 /||v (i, H 




= ("l) 1 1 p i llv (l, || 

i=l, . 

. . « k 

(9.8) 

x (i) = x (i) 

i=0, . 

. . , k 


Since the sets of vectors {v^ v^J and fr^, . . . , "r ^ are each 

mutually orthogonal it follows that the sets [v^, . . . , v^} and [u^ f i.., u^} 
defined by Equations (9.. 1) and (9.2) are each mutually orthonormal. Using 
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the equations of Algorithm CGC (5. 7) it can be directly verified that the 
quantities defined by Equations (9. 1) ~ (9.8) satisfy the relations of Algorithm 
ITC (8.31). Thus the algorithms CGC and ITC theoretically produce the same 
sequence of approximate solution vectors. 

Since the difference between the two algorithms only involves different scale 
factors it is to be expected that, apart from questions of exponent overflow or 
underflow, the two algorithms will exhibit essentially the same numerical 
behaviour also. 

The theoretical equivalence of these two algorithms was pointed out by 
Paige (1972), p. 13. 
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Chapter IQ Solving a Consistent System Ax=b where A is Symmetric 


Let A be an n x n symmetric matrix and let b be an n-vector contained in 
the column space (range space) of A. We wish to find a vector x satisfying 

(10. 1) Ax = b 


In practical problems of this type the matrix A would usually be of rank 
n. The algorithm to be described does not require that Rank (A) = n. If 
Rank (A) < n the solution of Problem (10. 1) is nonunique. In this case the 
algorithm finds the (unique) solution vector x lying in the row space of A. This 
is the minimal length solution vector for Problem (10. 1). 

Assume the existence of matrices V , and C, , and a k*-vector p such 

that 

(10.2) V = [v {1) ,...,v <k) ] 


C = 



(10.4) V T V = I k 

(10.5) AV = VC 


and 
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( 10 . 6 ) 


A „ - - A 

x = Vp 

If the matrices V and C are available Problem (10. 1) can be attacked 
as follows. Make the change of variables 

(10.7) x = Vp 

in Equation (10. 1), then use Equation (10. 5) obtaining 

(10.8) VCp = b 

T 

Left multiply Equation (10.8) by V using Equation (10.4). 

(10.9) Cp = V T b = g 


T 

Thus Problem (10. 1) could be solved by first computing g = V b, next 
solving Cp = g for p, and finally computing x - Vp. 

Since C is a symmetric tridiagonal matrix the entire matrix C must be 
determined before any components of p can be computed. Thus the equation 
Cp “ g is not directly suitable for use in an algorithm which discards old 
vectors v^ as new ones are computed. 

Let Q be a k x k orthogonal matrix which on post- multiplication times 
C produces a k x k lower tridiagonal matrix L. 
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'll 


°Zl *22 


( 10 . 10 ) 


CQ ~ L = 


b 31 


& 32 33 


0 £ k, k-2 ^k, k" 1 £ kk 


Define 


(10. 11) 


W = VQ 


and 


( 10 . 12 ) 


z = q t p 


Then 


(10. 13) 


Lz = CQQ T p = Cp = g 


and 


(10. 14) 


x = Vp » VQQ T p = Wz 


Thus Problem. (10. 1) can be solved by the following sequence of operations 
assuming that b } V, and C are known a priori 
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(10. 15) 

T 

g = V L b 

(10. 16) 

L- = CQ 

(10. 17) 

Solve L* = g for 2 

(10. 18) 

W = VQ 

( 10 . 19) 

x = Wz 

Each individual step of this sequence is amenable to being implemented in a 
sequential manner so that in fact V and C need not be known a priori but rather 
the column vectors of V and the elements o' . and 8^ of C can be computed, used, 
and discarded sequentially. 

This method of solving Problem (10. 1) is due to C, C. Paige and M. A. 
Saunders (Personal correspondence, 1972)„ 

We follow M. Saunders in defining and v^ by 

(10. 20) 

-3D 

►— * 

II 

and 


(10.21) 

v^ = b /p 1 [assuming p^O] 

From Equation (10. 5) we may write the equations 
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( 10 . 22 ) 


v^VAv^-v* 1 ^ 


(10.23) v (l+I) B i+1 = Av^ - v (l) *.- v* 1 " 1 ^. i=2, , . . , k- 1 


and 


(10.24) 


0 = Av 


(k) 


(k) 

- v' 'ey, 


(k ■ 
- v 


D c 


The numbers may be computed as normalization factors to assure 

/ i*f* 1 ) 

that the vectors v have unit length as required by Equation (10.4). The 
numbers c. can be computed as 

(10.25) o. = v (l)T Av (l) 

which is the condition (see Equation (2. 1)) which assures that computed 

by Equation (10.22) or (10.23) will be orthogonal to v^. 

The integer k will generally not be known in advance and may (theoretically) 
be determined as the first value of i for which the right side of Equation (10,23) 
is zero. With k so determined it can be verified that the vectors v^, . . . , v^ 
produced by use of Equations (10.21) - (10.23) form an orthonormal basis for 
a subspace of the column space of A. This verification makes essential use of 
the symmetry of A. 

It is pointed out by Paige and Saunders that the computation of the 
orthogonal set of vectors v^, i=l ,..., k, using Equations (10. 22)— (10. 23) is 
a method due to Lanczos (1950 and 1952). 
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The orthogonal matrix Q will be (implicitly) constructed as the product 


(10.26) Q = GjGg* • •G k _ 1 

whe re 


(10.27) G . = 


i- 1 2 k-i-1 


I 0 0 

0 G. 0 

l 

0 0 I 


i-1 

2 

k-i-1 


i— !»..-» k— 1 


(10. 28) 



i=l, . . . , k-1 


and 


(10.29) 


2 

c. 

l 



1 


Each matrix G^ effects a nontrivial transformation on a particular 3x2 
or 2 x 2 submatrix of the appropriate intermediate matrix which arises in the 
process of transforming the symmetric tridiagonal matrix C to the lower 
tridiagonal matrix JL* This action may be expressed as follows: 

(10. 30) T u = 0i l 
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(10.32) 


i= 1 1 . . . , k-2 


l. . 
1 

0 i+l 



0 

*i+i,i 

"i+l 

• G i = 

1 , i 

A i+1, i+1 

0 

B i+2_ 


_^i+2, i 

£ i+2,i+l_ 



— 


** — 

*k-l, k-1 

B k 


L , . i 0 

k" 1 1 k" i 

k-1 

V 

• G k-1 = 

_^k, k- 1 ^k, k_ 


As each additional row of L is determined by these equations an additional 
component z^ of the vector z satisfying Equation (10. 17) can be determined. Note 
that 


(10. 34) 


g - 


due to Equations (10.4), (10. 15), and (10.21). Thus the equations for the 
components z. are 


(10. 35) 


z i * Vhi 
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(10. 36) 


z 


2 = 


* 21*1 



and 


(10. 37) 


z. = -( i. . 

l i, i- 


,z. 

2 l 


■2 i, i-l 


i-l )/A ii i=3. 


k 


Define the vectors 


(10. 38) 



z i’ °> ■ ■ ' » » i-0, 1, . . . , k 

k-i~ 


and, motivated by Equation (10. 19), define a sequence of approximate solution 
vectors 


x <l} = Wz (l) i=0, 1, . . . , k 


The associated residual vectors are expressible as 


(10.39) 


r^ = b-Ax^ = b-AWz^ 


= b-AVQz*^ = b-VCQz* 1 * 


= b-VLz (l) = v^B^VLz^ 


. = VCe^^Bj-Lz^h 

from which we may write 


i=0, 1, . . . , k 
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(10. 40) 




(10.41) rt 1 ) = -[ v ( i+1 >. v (i+2) ] ‘ 


^i+1, i-l ^i+1, i 


vT 

0 . 

i+2, i_ 


— i 

•*4 


i=l, . . . » k-1 


(10.42) 


r< k > = 0 


Interpretation, of Equation (10.41) for i- 1 requires the definitions q = 0 anC ^ 
z Q =0 while for i=k~l one must define v ^ + ^=0 and = 

The norms of the residual vectors are expressible as 


(10.43) 

(10.44) 


P 0 s |M°*|| = 


p. s !lr^|i = [(X. , , . ,z. ,+X., , .z.) 2 +(X. , ,z.) 2 ] 

K i ' 1 L ' i+l, i-l i-l i+l,i i' ' i+2,i v 


1/2 


for i— 1 , . . . , k— 1 


(10. 45) 


0 k - l|r (k) ll = o 


where again we define q=0, Zq=0, and Ic-1 = ®‘ 

Combining these equations appropriately one can obtain the following 
algorithm 
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(10.46) Algorithm ICSE Iterative Solution of a Consistent Symmetric System 

of Linear Equations [Originated by C. C. Paige and 
M. A. Saunders, personal communication, 1972^ 


Step 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 


Description 

*<°W 

S,: = IM 

If Bj-0 set k: = 0 and go to Step 24 

v^^b/pj 

u (U :=v (D 

QUIT:= FADSE 
i: = l 

y (i Wv (i) 

(i)T (i) 
a.:=v y 


~(i+l). 
v : 


y* 1 *-*/'* ifi=l 
y^-Of.v^-0.v^ ^ if i > 1 


11 B i+1 H!v (i+I >ll 

12 Theoretical termination test: If 3 i+ ^ = 0 set QUIT := TRUE 
Practical termination test: If B. + 1 is sufficiently small 


set QUIT: = TRUE 



i 


7 11 

It 

*1 

<ni 

1 


_ 6 2 _ 


4. . , 
1 , i-I 

j e. . 


l. . | O'. 

1, 1- 1 1 



1 

M* 
+ 
*— ■ 

H‘ 

1 

H- 

x., . 

i + i , i_ 


_° h + l 


v 1 - 
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Step 


14 

15 

16 

17 

18 


19 

20 

21 

22 

23 

24 


Description 

/ 

If i- 1 go to Step 1 6 

/ If i=2 set pji * -(^ 21 z 1 ) 2 +(^ 31 z 1 ) 2ll/2 

If i >2 set p. . 0 z, + £, , z. ,) 2 +(f.,, . ,z. 

^ K i-1 1,1-2 i-2 i,i-l i-l' ' i+I, i-I i- 

Practical termination test: If ^ is sufficiently small 

set k: —i-l and go to Step 24. 

— 2 2 1/2 
+ B tV 


^ 1 / ^ 1 1 ^ i "* 


Z i : " i ”' e 21 z 1^22 *~ 2 

-(JL. . 0 z. j+l. . ,z. ,)/&.. if i > 2 

^ ' i, i-2 i-2 i, l- 1 i- 1' ii 


If QUTT= FALSE set v (l+1) :=v {l+1 ] /B. ,, c.: = i../i.., 

i+l l ii ii 


V^iH^ii’ 



andrw (i) ,u (i ‘ fl) l: = [u (i) ,v (i+1, ] G . 


If QUIT^ TRUE set 

x (i W i_1) +w (i) z. 

l 

If QUIT= TRUE set k:=i and go to Step 24 

'-‘t " 

i: = i+l 

Go to Step 8 

(k) 

Here the algorithm is finished with the solution vector x' . 
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